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Abstract 

The process of genetic recombination can be seen as a chemical reaction 
network with mass-action kinetics. We review the known results on 
existence, uniqueness, and global stability of an equilibrium in every 
compatibility class and for all rate constants, from both the population 
genetics and the reaction networks point of view. 


1 Introduction 


Recombination, or chromosomal crossover, is the exchange of genetic mate¬ 
rial between a pair of homologous chromosomes. It occurs when matching 
regions on matching chromosomes break and then reconnect to the other 
chromosome. Recombination is the major evolutionary force to produce 
and maintain variation in a (sexual) population. 

The process of recombination 
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involves two loci, two alleles per locus, and hence four gametes. It can be 
written as a reversible chemical reaction 

9i + 9 a & 92 + 93- 
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This is a very simple reaction network with deficiency zero. The recom¬ 
bination 51+54 —>■ 52 + 53 occurs at the rate kpipA determined by the 
rate constant times the frequencies of the reacting gametes. In the chemical 
setting, this corresponds to the assumption of mass-action kinetics. The 
dynamical system for the gamete frequencies amounts to 

Pi=P 4 = k (~PlP 4 + P2P3) = ~P 2 = P 3 - 

There are conservation laws 


Oi +P2)' = {pi +P3)' = {P1+P2)' = {P1+P3)' = 0 

of which three are linearly independent. The equilibrium manifold is given 
by the conic 

PlPA = P 2 P 3 - 

In each stoichiometric compatibility class, solutions converge to a unique 
detailed-balancing equilibrium. This is usually proved by considering the 
so-called linkage disequilibrium function D = P1P4 — P2P3, which satisfies 

D = -kD (pi +P2+P3+ Pa) ■ 


Hence, D converges to 0. 

Alternatively, one can try an ansatz for a Lyapunov function in the form 

V{p) = Y,Li p {Pi)- Then, 

4 

V(p) = J2 F '{Pi)Pi = ~kD(F\ Pl ) + F\ P a) - F\p2) - F 1 (p 3 )) . 

i=l 

For the choice F'(p) = In p, we obtain 

V(p) = —kD (Inpi + lnp4 — lnp2 — lnp 3 ) 

= -k (piPA - P2P3) (ln(piP4) - ln(p 2 p 3 )) < 0 

due to the monotonicity of the logarithm. This shows that V(p), with the 
convex function F(p) = plnp — p, is a global Lyapunov function. 

The main object of this paper is to study the general recombination model in 
continuous time, with an arbitrary number of genetic loci and arbitrary num¬ 
bers of alleles at each locus. We will see that this leads to a chemical reaction 


network which is reversible, satisfies the Wegscheider conditions 22j] (since 
the rate constants of a reaction and its reverse coincide), and is detailed- 
balancing. A generalization of the above entropy-like Lyapunov function 
allows us to prove global stability. 

In population genetics, this general recombination model was studied (in 
discrete time) by Geiringer 0, and further by Their proofs use 

linkage disequilibrium functions, induction on the number of loci, cumulants, 
etc., and are far from easy. Simpler proofs based on the entropy as Lyapunov 
function were independently given by Akin fl] (in continuous time) and 
Lyubich [2| (in discrete time). 
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2 Mathematical model 

Notation: We denote the positive real numbers by R>o and the non¬ 
negative real numbers by R>o- For a finite index set I, we write M 1 for 
the real vector space of formal sums x = Yliei x i * with Xi £ R. Viewing 
the elements of I as indicator functions, a vector x £ R^ can be seen as a 
function x: I — > R, and x(i) = X{. For x,y £ R> 0 , we dehne x y £ R>o as 
where we set 0° = 1. Given a matrix Y £ R /xJ , we denote by 
W £ R 7 the column vector indexed by j £ J. For x £ and Y £ R>q J , 

we define x y £ R> 0 as (. x Y )j = x YJ = n, e / x J' :i for j £ J. 


2.1 Genetic recombination 


We consider a finite set of loci C with L = \C\ > 1 , finite sets of alleles A% 
with Ai = \Ai\ > 2 for i £ £, the resulting set of gametes 

fy = x ... x dh, 

and the set of recombination patterns 

V = {{/, J} | / c C, J = C\I}. 

In a recombination following pattern {/, J}, alleles at loci I get separated 
from alleles at loci J. There are \Q\ = Ai gametes and \V\ = 2 L_1 

recombination patterns, including the trivial recombination { 0 , >C}. 

Further, we introduce the distribution of gamete frequencies 

P- G -» K > 0 

and the distribution of recombination rate constants 


c.V^f M>o- 

We identify the function p: Q —>• R>o with the vector p £ R > 0 and write 
p = YlgegPid) 9■ Usually, we are interested in elements of the simplex 

Sg = {p£ M | 0 | Y. g& gP{9) = !}■ 


For gametes g,h £ Q and a recombination pattern {/, J} £ V, we dehne 
gihj £ £/ as 


G 9ihj)i 


9i, Ai £ I, 

hi, if i £ J 


and the resulting recombination as 


{ 5 , h} 






( 1 ) 
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with rate constant c({/, J}). For gihj, gjhj £ Q and {I, J} £ V , we find 
(gihjjiigjh^j = g and (gihj) = h 

and obtain the reverse recombination 

{gihj^jhj} c({ i/ }) {^/i} 
which occurs with the same rate constant. 

Clearly, recombination © causes a change in gamete frequencies only if 
{g,h} 7 ^ {gihj, gjhj} and c({/, J}) > 0. Further, different recombination 
patterns may give rise to the same recombination (with different rate con¬ 
stants, in general). In order to view recombination as a chemical reaction, 
we have to ensure inequality of left- and right-hand sides and positivity of 
rate constants. Moreover, we have to sum over the rate constants of all con¬ 
tributing recombination patterns which can be seen as reaction mechanisms. 

2.2 Chemical reactions 

Let K C C be a subset of loci. The recombination pattern {/, J} £ V 
induces the subpattern {I, J}k £ Vk where {/, J}k = {I D K, J fl K} and 

V K = {{/, J} | I C K, J = K \ I}. 

We write {/, J} > {I,J}k and, for simplicity, V* K = Vk \ {{0,iF}}. The 
set of all recombination subpatterns amounts to 

V = (J V K , 

KCC 

and we introduce the distribution of cumulative recombination rate con¬ 
stants 


C'.V^r M>o, 

{} 

For a subpattern {/, J} with J = K \ I and K C C , the cumulative rate 
constant c({/, J}) sums over all patterns which agree with {/, J} on K. 

An important parameter is the cumulative rate constant for the recombina¬ 
tion subpattern {{i}, {j}}, that is, for the case that an allele at locus i gets 
separated from an allele at locus j. We assume that c({{i}, {j}}) > 0 for 
all pairs of loci i,j. Otherwise, the two loci can be identified. 
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To explicitly state a chemical reaction arising from a recombination pattern 
and a pair of gametes, we introduce the set A (g, h) = {i £ C \ gi ^ hi} for 
gametes g,h £ Q. In genetic terms, g and h are heterozygous at the subset 
of loci A(g,h) and homozygous otherwise. 

Now, gametes g,h £ Q and a recombination pattern {/,«/} £ V give rise 
to a reaction mechanism, only if \A(g,h)\ > 2, {I, J}A( g ,h) / {0, A(g, h)}, 
and c({I, J}) > 0. In other words, only if the gametes are heterozygous at 
two or more loci, if the subpattern induced on these loci is non-trivial, and 
if the recombination rate constant is non-zero. Every pattern {I', J’} £ V 
with {/', J'} > {I, </}a( 3 ,/)) and c({I', J'}) > 0 gives rise to a mechanism for 
the same reaction, that is, to the same recombination (with different rate 
constant, in general). The effect of all such patterns is summarized in the 
chemical reaction 

k 

9 + h > gihj + gjhj (2) 

with rate constant k = k(g + h —> gihj + gjhj) = c({/, J}A(g,h)) > 0. Note 
that g + h stands for {g, h} such that g + h equals h + g. The reverse reaction 

k 

9ihj + gjhj g + h 
occurs with the same rate constant. 

2.3 Reaction networks 

A chemical reaction network ( S,C,1Z ) consists of three finite sets: a set S 
of species, a set C C Mf 0 of complexes, and a set 1Z C C x C of reactions. 
Complexes are the left- and right-hand sides of reactions. A complex y £ C 
can be seen as a formal sum of species y = Ylsesys s i where y s is the 
stoichiometric coefficient of species s. For a reaction (y, y') G 1Z, we write 
y -+ y' • It is required that each complex appears in at least one reaction 
and that there are no reactions of the form y —> y. 

A chemical reaction network ( S,C,1Z ) together with a vector of rate con¬ 
stants k £ M> 0 gives rise to a weighted directed graph with complexes as 
nodes, reactions as edges, and rate constants as labels. The connected com¬ 
ponents of this graph are called linkage classes. (Note that linkage classes 
have nothing to do with genetic linkage.) A network is called weakly re¬ 
versible if every component is strongly connected, that is, if there exists a 
directed path from each node to every other node in the component. 

In the process of genetic recombination, the reacting species are the gametes, 
that is, S = Q. Every complex g + h is a formal sum (with stoichiometric 
coefficients equal to one) of two gametes g and h, which differ at two or more 
loci, and every reaction g + h —> gihj + gjhi arises from a pair of gametes 
and a recombination pattern {/, J}, under the conditions specified in the 
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previous subsection. The set of all chemical reactions (with corresponding 
rate constants) amounts to 


1Z = 


\g + h-+ gihj+gjhj 


g, h G Q, {/, J} G V with | A (< 7 , h) | > 2, 


{I,J}A( g ,h) ^ {0, A(g,h)}, and 
k = c{{I,J} HgM ) > o|. 


(3) 


For each reaction y ^ y' we have its reverse y 1 —> y, and both occur with 
the same rate constant. Hence, we can combine them in the reversible 
reaction y y', which we identify with y' y, and write k(y y') for 
k(y —> y') = k(y' —> y). From (HD, we obtain the set of all reversible reactions 


lZ^ = {y&y'\(y A y') € 1Z} 


(4) 


and the set of all complexes 


c = {y I (y y') e iz}. 


(5) 


In the examples and schemes below, we determine the set of all (reversible) 
reactions in another way. We first iterate over subsets of two or more loci 
and then over non-trivial subpatterns on these loci: 


7Z^> = 15 + ^-^ gihj + gjhj K C C, with \K\ > 2, 

g,h € Q with |A (g, h) \ = K, 

{I , J} e V* AM with k = c({/, J }) > o}. 


Thereby, we extend the definition of gihj to the subpattern {I, J} G V^g^) 
in the obvious way: (gihj)i = (ji for i G I, (gihj)i = hi for i G J, and 
= gi = hi for i G C \ (I U J). 

Finally, we consider the graph arising from the reaction network ( S,C,1Z ), 
in particular, its linkage classes. We observe that species (gametes) consist 
of alleles and complexes (pairs of gametes) contain two alleles at each lo¬ 
cus. Since reactions separate alleles, but do not consume or produce them, 
only complexes which contain the same alleles are connected by a reaction. 
Moreover, if complexes g + h and g' + h' are connected by a reaction then 
A(g,h) = A(</,/i'), and every subpattern {I, J} G V*^ gh ^ which gives rise 
to a reaction involving g + h gives rise to a reaction involving g' + h', and 
vice versa. Hence every linkage class is a symmetric graph. If no reaction 
is precluded by a zero rate constant, then every linkage class is a complete 
graph, characterized by two (possibly identical) alleles at each locus. 
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2.4 Examples and schemes 

We consider examples of genetic recombination for small numbers of loci and 
alleles and depict the corresponding chemical reaction networks as graphs. 
Further, we present schemes for arbitrary numbers of loci and compute the 
resulting numbers of linkage classes, complexes, and reversible reactions. For 
simplicity, we assume that no reaction is precluded by a zero rate constant. 
In this case, all linkage classes are complete graphs. 

Instead of c({J, </}) we write c(I) and further omit the set brackets, e.g., 

C ({{1},£ \ {!}}) = c ({ 1}) = c (l)- 

Example 1 (L = 2 loci with A\ = A 2 = 2 alleles). 

0 , 1 c(i) 1 0 

+ < > + 

0 1 0 1 

The graph has l = 1 linkage class, m = 2 complexes, and r = 1 reversible 
reaction. 

Example 2 (L = 3 loci with A\ = A -2 = A 3 = 2 alleles). 
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c(l)+c(2) 




1 0 


0 1 

0 + 0 

5 

1 + 1 

0 1 


0 1 

0 0 


1 1 

1 + 0 

5 

0 + 1 

0 1 


0 1 

0 1 


1 0 


c(l) 


1 + 0 


1 + 0 

0 1 

c(3) 

0 1 



0 


1 

c(3) 

1 


0 

0 

+ 

1 

c(l/ 

0 

+ 

1 

0 


1 

0 


1 



The graph has l = 7 linkage classes, m = 16 complexes, and r = 12 reversible 
reactions. The last class has 2 i_1 = 4 complexes and ( 2 2 ) = (.)) = 6 
reactions. 
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Scheme 3 (L > 2 loci with Aj = 2 alleles, i = 1, ..., L). 


< = £ 


->A-fc 


k=2 


= 3 L - (2 l + L 2 L ~ l ) 
= 3 l — 2 l ~ 1 (2 + L) 


L 

m = Y.] 

k =2 



2 -£/—k 2&-1 


= 2 i " 1 ( 2 i - (1 + L)) 



= EL 2 (£) 2 i_fc 2 fc_1 (2 fc_1 — 1) 2 _1 

= EL 2 (fe) (2 L_3 2 fc — 2 l ~ 2 ) 

= 2 l ~ 3 (3 l - (1 + L 2)) - 2 l ~ 2 (2 l - (1 + L)) 

= 2 L " 3 (3 i - 1) - 2 L " 2 (2 i - 1) 

Example 4 (L = 2 loci with Ai = 2 and ^2 = 3 alleles, respectively). 
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1 

c(l) 
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0 

0 

+ 

ll 

' X x * 

0 

+ 

1 
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1 

c(l) 

1 


0 

0 
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X- 

0 
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2 


2 

0 


1 

C(l) 
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0 

1 

+ 



1 

+ 


2 


2 


The graph has l = 3 linkage classes, m = 6 complexes, and r 
reactions. 

Scheme 5 (L > 2 loci with .Aj > 2 alleles, i = 1, ..., L). 


<- e n 

irC£:|A'|> 2 ieAT 

m= e n 

A'C£:|A'|> 2 ieAT 

- e n 

A'C£:|A'|> 2 i£AT 


n ^ 

iec\K 

n A i2 |Ar|_i 

i&C\K 

n a 

iec\K 


3 reversible 
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2.5 Dynamics 

Recombination P) causes a change in gamete frequencies proportional to 
gihj + gjhi — g — h at the rate c({7, J}) p(g) p(h) determined by the re¬ 
combination rate constant times the frequencies of the “reacting” gametes. 
We formulate the dynamical system for the vector p £ of all gamete 
frequencies, that is, for p = J2 g eg P(d) 9i by summing over all recombination 
partners and patterns: 

= \ c({I,J})p(g)p{h)(g I hj + gjh I -g-h). ( 6 ) 

g,h&g {i,j}ev 

The contribution of a particular recombination m is identically zero if 
{g, h} = {gihj, gjh]} or c({/, J}) = 0. On the other hand, different re¬ 
combination patterns may yield the same recombination (except for the rate 
constants). As detailed in Subsection 12.21 the effect of all patterns causing 
recombination P) can be summarized in the chemical reaction P, provided 
that the recombination is effective and the cumulative rate constant is pos¬ 
itive. Reaction @ occurs at the rate c({I, J}^ g ^)p(g)p(h). In chemical 
terms, it follows mass-action-kinetics. 


Mass-action kinetics 

Let (S, C, TV) be a chemical reaction network and k £ R> 0 a vector of rate 
constants. Under the assumption of mass-action kinetics, the rate of a reac¬ 
tion (y —»• y') £ 77. which depends on the species concentrations x £ M> 0 , is 
given by k{y —> y') x y , that is, by a monomial in the reactant concentrations 
with the corresponding stoichiometric coefficients as exponents. 

Hence, we obtain a dynamical system equivalent to ([ 6 ]), by summing over 
all reactions ([3]) and assuming mass-action kinetics: 

% + h ->• g' + ti) p(g) p(h ) (g 1 + h! - g - h) . (7) 

(g+h—tg'+h')£lZ 

The right-hand side of © can be written as a product of the stoichiometric 
matrix N £ and the rate vector Vf.(p) £ Thereby, the column 

vector of N indexed by (g+h —> g'+h') £ 77 is given by (g' + h' — g — h) £ 
and the component of Vk(p) indexed by g + h —> g' + h! is given by k(g + h —»• 
g 1 + h!)p(g)p(h). Hence, 


dp 

dt 


Nv k {p). 


( 8 ) 
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Complex balancing 

The right-hand side of the dynamical system 0 can also be written as a 
product of the complex matrix Y £ M^ xC , the Laplacian matrix £ M CxC 
of the weighted directed graph, and the vector of monomials p Y £ The 
column vector of Y indexed by y £ C is given by y € Rg 0 itself, that is, 
Y y = y, and Ak is defined as follows: ( Ak) y ' y = k y ^ y ’ if (y —>• y') £ 1Z, 
( A k ) yy = -J2( y ^y')e-R. k y^-y'’ and ( A k)y'y = 0 otherwise. We obtain 

Recall that (p y ) y = p' 1 v = p y for y £ C. For a particular complex y = g + h, 
we have p y = p 9+h = p(g)p(h). 

An equilibrium of (j9]) is called complex-balancing if A^p 1 = 0. That is, if 
at each complex the rates of all reactions sum up to zero. 

Detailed balancing 

In the process of genetic recombination, all reactions are reversible. More¬ 
over, the rate constants of a reaction and its reverse coincide. Hence, we 
obtain a dynamical system equivalent to CD, by summing over all reversible 
reactions ( 0 : 

^ K 9 + h&g' + ti)(p(g)p{h)-p(g')p{ti)){g' Ah' -g-h). 

(g-\-hOg' -\-h')£ TZo 

(10) 

An equilibrium of (HOl) is called detailed-balancing if p{g)p(h) = p{g')p{h') 
for all (g + h g 1 + h') £ IZ^,. In general, an equilibrium of a reversible 
reaction network is called detailed-balancing if the rates of each reaction 
and its reverse coincide. Clearly, every detailed-balancing equilibrium is 
complex-balancing. 

2.6 Conserved quantities 

The change over time 0 lies in a subspace of M®, and every trajectory in M > 0 
lies in a coset of this subspace. We define the stoichiometric subspace 

S = irn A r 

and the stoichiometric compatibility classes 

S( P ) = (p + S)nR g > 0 

for p £ M>q. Every stoichiometric class is characterized by its orthogonal 
projection on S Y = (imIV)- 1 - = kerIV T , that is, by a vector of conserved 
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quantities. For u G S'" 1 , that is, u T N = 0, we have 

d{u T p) 
d t 

that is, u T p = const. 

We observe that the vector 1 = I s = J2 g eg 9 is orthogonal to all columns 
of N: 1 T (g' + b! — g — h) = 0 for all (g + h — > g' + h') € 1Z, that is, 1 T N = 0. 
Since 1 T p = J2 g£ gP(9), we have 

d(E 3 egP(5)) 

d t ~ ’ 

and, as one consequence, the simplex Sg is invariant. 

Further, we consider for each locus and each allele at this locus the subset of 
gametes which contain this allele and define the corresponding formal sum 
of gametes 

Ui(a ) = g for i G C and a G A x . 

geg.gi=a 

where iq(a) G {0,1}^. As already mentioned, only complexes which contain 
the same alleles are connected by a reaction. Hence, Ui(a) T (g'+h'—g—h) = 0 
for all (g + h —> g' + h!) G 1Z, that is, Ui(a) T N = 0, and the marginal 
frequencies 

Pi{a) = Ui{a) T p = ^2 Pig) 

g&Q-gi=a 

are conserved quantities, that is, 

dPi(o) = Q 
d t 

For each i G £, we have 

Ui (°) = Yl , 9 

and 

J2 ^( a ) = ^2p(p)- 

Hence, there are at least l + ^A e £(Aj —1) linearly independent vectors in S L 
and as many independent marginals. 

We define the marginal compatibility classes 

M(p) = {p G M> 0 | Pi(a) = Pi(a) for i G £ and a G Ai} 
for p G M> 0 . Clearly, S(p) C M(p). 
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3 Results 

We determine the equilibria for the process of genetic recombination and 
prove convergence to a unique equilibrium. 

First, we rewrite the dynamical system ([6|) . Using the symmetry in the 
double sum over recombination partners, we obtain 

= E E C (R J })p(g)p( h ) idihj-g) 

g,heQ {i,j}ev 

= E E p(g)p(h ) gihj-p 

{i,J}eP \g,heg 

Thereby, we assumed E 9 g gP(g) = 1, that is, p £ Sg. 

For K C C, we define the set of subgametes Qk = FLeA'^*’ the projection 
Q Qk, 9 *->• gK, where {gx)i = 9% for i £ K, and its linear extension 
to the corresponding vector spaces: R> 0 —> P *->■ Pk, where pk = 

ZgesP(9)9K, that is, p K {9K ) = E heg-.h K =g K pW- If K = I 2 '} with i G C, 
we recover the marginal frequencies pi = E g ^gP(g)9i, that is, Pi(gi) = 

E h&g-.hi= gi PQ 1 ). 

Let {I, J} £ V. For g £ Qi and h £ Qj, we define gh £ Q as (gh), = gi for 
i £ I and (gh), = h, for i £ J and extend the multiplication Qj x Qj^Q 
linearly to M> 7 0 x —> M> 0 . Hence, we write 

E C ({ J ’ J }) \^2p(g)gi ^2p( h ) hj -p ] ( n ) 

{I,j}er \geg h&g ) 

= E c ({ J ’«/}) (pipj -p) ■ 

{i,J}eV 

In fact, we may sum over {I, J} £ V* since the contribution of {0,£} is 
identically zero. 

The projection of a trajectory of the dynamical system is the trajectory of 
a projected dynamical system with the same structure: For K C C and 
{I', J'} = {/, J}k G Vk, we find 

(pipj)k = E p(a)p(g) (gi h j)K = E p(g)p(g)9i ,h J' =pvpj> 

g,heg g,h&g 
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and hence 

= X X C ( U ’ J }) ((pipj)k ~Pk) 

{i',J'}ev K {i,J}eP: 

= X c({l', J'}) ( p P pj' - p K ) ■ 

3.1 Equilibria 

Now, we are in a position to characterize the equilibria on the simplex. 

Lemma 1. For all recombination rate constants, p € Sg is an equilibrium 
of the dynamical system ((6j) if and only if 

p = lift, that is ’ p (s) = II«(&)■ ( 12 ) 

iec iec 

Proof. We show that, if (fT2|). then 

PI = W Pi 

iei 

Indeed, for {/, J} € V and J = {1,..., | J|}, we find 

P/(S/) = X] P^ lh j) 

hjSGj 

= x ■■■ x nft<*)iift(fc) 

hj.se! h\j\SQ\j\ isl iSJ 

=n«(ft) x] ft^ 1 ) ••• x] pijk^iji) 

zG/ /i|j|G^|j| 

=nftCft)- 

iSl 

Hence, 

PiPJ = lift lift = lift =p 

iSl iSJ iSC 

for all {/, J} G V, and p £ Sg is an equilibrium of the dynamical system (Hill 
equivalent to ©. 

It remains to show that p £ Sg is an equilibrium only if (1121) . We proceed 
by induction on the number of loci: 

For C = {1}, there is no non-trivial recombination. Every p € Sg is an 
equilibrium which coincides with its marginals: p(g) = pi(g) for g £ Q, that 
is, p = p\. 
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For L > 2, we consider subsets of loci K C C with \K\ < L. The projection 
of an equilibrium p £ Sg of the dynamical system (with loci C) is an equi¬ 
librium of the projected dynamical system (with loci K). By the induction 
hypothesis, px = ThexPi an d hence 

PiPJ = \\pi Pi = in* 
iei i&J iec 

for all {I, J} £ V*. Summing over {I, J} £ V* in (fill) , we obtain 

0= c ({^ J })\YlPi-P 

{i,J}&V* \ieC 

and hence p = Y\i^gPi- (There always exists a non-zero rate constant for 
some non-trivial recombination pattern.) □ 

If the dynamics is not restricted to the simplex, then p £ M> 0 is an equilib¬ 
rium if and only if 

p= (EgegPte)) lift 

ie£ 

Clearly, every equilibrium is determined by its marginals, and we have the 
following result. 

Proposition 2. Every marginal compatibility class contains a unique equi¬ 
librium. 

In chemical terms, every equilibrium is detailed-balancing since 

p(g)p{h) = Y\pi{gi) Pi(hi) = J \pi(g'i) Pi(h'i) = p{g')p{h') 
iec iec 

for all (g + h <=> g' + h') £ 1Z&., cf. Equation (fTOT) . Recall that only complexes 
which contain the same alleles are connected by a reaction. In fact, we can 
derive the following result entirely in the chemical setting. 

Proposition 3. Every stoichiometric compatibility class contains a unique 
positive equilibrium, which is detailed-balancing, and no boundary equilibria. 

Proof. First, we determine the set of positive detailed-balancing equilibria. 
Using positivity, we write the condition for detailed balancing, 

P{g) P{h) = p{g') p(h') for all (g + h g' + h') £ K^, 

as 

po'+h'-g-h = x for all ( g + h -£ g ! + ti) £ U 
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and even more abstractly as 

= i, 

where N G is the stoichiometric matrix and 1 = l 7 '’. Clearly, the 

trivial solution is given by p* = 1 = 1^. To determine all solutions, we take 
the logarithm, 

N t In p = 0, 

and note that ker N T = (imJV) 1 = S 1 - and dim(S'" L ) > 1. Hence, we can 
write the set of positive detailed-balancing equilibria as 

Z = {p £ M> 0 | In p — lnp* G 5 X }. 


Clearly, every detailed-balancing equilibrium is complex-balancing. Now, 
if there exists a positive complex-balancing equilibrium p*, then the set of 
positive complex-balancing equilibria is given by Z and there are no other 
positive equilibria 1CJ, Theorem 6A], Moreover, every stoichiometric compat¬ 
ibility class contains a unique positive equilibrium [10, Lemma 4B], Hence, 
the sets of positive detailed- and complex-balancing equilibria coincide, and 
every stoichiometric compatibility class contains a unique positive equilib¬ 
rium, which is detailed-balancing. 


It remains to preclude boundary equilibria. We consider an arbitrary initial 
value p G K|„ on the boundary, that is, p(g) = 0 for some g G G, and 
define the set of gametes Qq = {g G Q \ p{g) > 0}. Note that for each locus 
and each allele at this locus there is a gamete which contains this allele and 
occurs with a positive frequency. Further, recall that each pair of loci gets 
separated by some recombination pattern with positive rate constant. By 
Lemma [H below, the set of all gametes Q is reachable from Qq. Now, let p(t) 
be the solution of the dynamical system with p(0) = p. By [2l|, Theorem 2, 
p. 618], p(t) G M| 0 for t > 0, and hence p is not an equilibrium. □ 


In the proof, we have used fundamental results about complex balancing 
by Horn and Jackson [10]. Necessary and sufficient conditions for complex 
balancin g ar e given by Horn [9 ! ] and for detailed balancing by Vol’pert and 
Hudjaev 21] and Feinberg 0]. For the relation between complex and detailed 
balancing, see 


In the proof of Proposition (3[ we have also used the purely graph-theoretical 
concept of reachability. Let S and 1Z be the species and reactions of a 
chemical reaction network, and let <S>o C S be a set of species. Iteratively, 
we define 


Si — Si —i U {g | g, h G Si -1 and (g + h -G g + h!) G 7Z} 

for i > 1. Since the graph is finite, we find S t = Si * for some i* > 0 and all 
i > i*, and the set of species reachable from <Sq is given by Si*. 
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For chemical reaction networks arising from genetic recombination, we have 
the following result. 

Lemma 4. In a process of recombination, assume that every pair of loci 
gets separated by some pattern with positive rate constant. In the resulting 
chemical reaction network, every gamete is reachable from a given set of 
gametes, provided that every allele is contained in some gamete in this set. 

Proof. We use induction on the number of loci: 

For L = 1, the gametes coincide with the alleles. 

For L > 2, we consider subsets of loci K C C with \K\ < L and project the 
network and the given set of gametes on the loci K. In the projected net¬ 
work, every pair of loci gets separated by some pattern, and in the projected 
set, every allele (at loci K ) is contained in some gamete. By the induction 
hypothesis, every gamete gx € Gx is reachable in the projected network, 
and hence some gamete h € G with hx = 9k is reachable. 

It remains to show that every gamete g' € G is reachable. Let C = {1,..., L} 
and G = £\{1}, H = £\{2}. By the argument above, some gametes g,h £ G 
with go = g' G , hn = g'x are reachable. If g' equals either g or h , then it is 
reachable. Otherwise, since loci 1 and 2 get separated by some pattern, we 
find the reaction g + h —>• g' + hi , and hence g' is reachable. □ 

On the one hand, by Proposition [2J every marginal compatibility class con¬ 
tains a unique equilibrium. Hence, the set of all equilibria can be parame¬ 
trized by 1 + — 1) independent marginals. On the other hand, 

by Proposition O every stoichiometric compatibility class contains a unique 
equilibrium. Since stoichiometric classes are contained in marginal classes, 
we have the following result. 

Corollary 5. The stoichiometric compatibility classes coincide with the 
marginal compatibility classes, and dim(S’- L ) = 1 + ~ !)• 

3.2 Convergence 

Our main results concern the convergence of the dynamics to a unique equi¬ 
librium. For the first theorem, we provide two proofs: one by induction (as 
in the original literature) and one using the entropy as a Lyapunov func¬ 
tion. For the second theorem, formulated in the chemical setting, we rely on 
results from chemical reaction network theory which are based on the same 
Lyapunov function. 

Theorem 6. In every marginal compatibility class and for all recombination 
rate constants, a process of genetic recombination converges to the unique 
equilibrium given by m- 
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First Proof (Induction). Every marginal compatibility class is characterized 
by a unique equilibrium. Given an equilibrium p £ Sg, that is, p = H iecPh 
we consider trajectories in the class M(p). We proceed by induction on the 
number of loci: 

For L = 1, we have M(p) = p. 

For L > 2, we consider subsets of loci K C £ with \K\ < L. The projection 
of a trajectory (f> : R>o — > M(p) of the dynamical system (with loci C) 
is a trajectory of the projected dynamical system (with loci K). By the 
induction hypothesis, <f(t)K —> II ieKPi as t —> oo and hence 

<t>(t)j -^Y[piY[pi = Y{pi=p 

i£l i£j iec 

for all {/, J} £ V*. Summing over {I, J} £ V* in the dynamical system (HID 
equivalent to ([6]), we obtain the non-autonomous differential equation 

T7= C ({ J ’ J}) (<l>i~ $) 

{i,J}ev* 


with 

f(t)= C ({ J > J })^^ 

and 

/(*)-»• 5Z 


as t —>• oo. In other words, the differential equation is asymptotically au¬ 
tonomous. The limiting equation 


{U}eV* 


is linear, and hence <f(t) —> p in the limiting equation. Moreover, {p} is 
the maximal compact invariant set in the limiting system, and therefore 
<j>(t) —>• p holds also for all solutions of the original dynamical system, see 


e.g. [13, 14]- 


□ 


Second Proof (Lyapunov function). We consider the classical entropy func¬ 
tion 

H{p) = -^2p(g)\np(g) = -p T lnp > 0 
g&G 
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which defines a continuous function on the simplex Sg. If p(g ) > 0 for all 
g £ G, then H is smooth and 

H(p) = - ^2p(g)\np(g) - ^2fi(g) = ~p T hip, 

g&G g&G 

since "Yh g& gp{g) = 1. Using the dynamical system (fTOl) equivalent to JB|, we 
obtain 

H{p)= k(g + h gf + h!) (p(g)p(h) -p{cf)p{h')) ■ 

(g+hog'+h ') 

• (In (p(g)p(h)) - In (p(g')p(h'))) > 0. 

Equality H(p) = 0 holds if and only if p is a detailed-balancing equilibrium, 
that is, if and only if (1121) holds. 

Given an initial point p(0) £ K?o in the interior, the entropy H(p(t)) in¬ 
creases strictly towards its maximum on M(p( 0)), and pit) converges to the 
unique equilibrium p in the class M(p( 0)). 

Given an initial point p(0) £ on the boundary, we have p(t) £ 
for t > 0, cf. the proof of Proposition [3j In genetic terms, recombination 
immediately produces all gametes, as long as all alleles are present in the 
population. □ 


The entropy as Lyapunov function was used by Akin [lj and Lyubich 12] 
(referring to a paper by Kun and Lyubich 0 ) to prove global stability for 
recombination. For chemical reaction networks with detailed balancing, see 
Vol’pert and Hudjaev 2(^2l| (who acknowledge previous use of the entropy 
function by Zel’dovich [231]). For complex balancing, see 0 S E0. 


Theorem 7. A process of genetic recombination gives rise to a reversible 
chemical reaction network with mass-action kinetics. In every stoichiomet¬ 
ric compatibility class and for all reaction rate constants, the dynamics con¬ 
verges to a unique positive detailed-balancing equilibrium. 


Proof. By [2ll . Theorem, pp. 642-643], the cu-limit set of every solution con¬ 
sists either of a unique positive detailed-balancing equilibrium or boundary 
detailed-balancing equilibria. By Proposition [3] there are no boundary equi¬ 
libria, and every solution converges to a unique positive detailed-balancing 
equilibrium. □ 


4 Final remarks 

Note that we have not used a central concept of chemical reaction network 
theory, the deficiency 


5 = m — l — s, 
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where m is the number of complexes, l the number of linkage classes, and s 
the dimension of the stoichiometric subspace. 

The deficiency zero and one theorems state that there exists a unique (asymp¬ 
totically stable) positive complex-balancing equilibrium, for all reaction rate 
constants and all stoichiometric compatibility classes, if the network is weakly 
reversible and either (0) its deficiency is zero or (la) the deficiencies of the 
individual linkage classes are zero or one and (lb) the individual deficiencies 
add up to the total deficiency, see [5]. 

In Example 1 (L = 2, A\ = A 2 = 2), we find <5 = 2 — 1 — 1 = 0. However, 
already in Example 2 (L = 3, A\ = A 2 = A 3 = 2), the deficiencies of the 
individual linkage classes are zero, but <5 = 16 — 7 — 4 = 5. 

In fact, the individual deficiencies are zero in the entire Scheme 3 (L > 2, 
Ai = 2 for i = 1,..., L): Every linkage class is characterized by two different 
alleles at some loci K € C with \K\ > 2 and two identical alleles at other 
loci. Hence 2l A l~ 1 — 1 — (2 ^ K ^ 1 — 1) = 0, whereas 


<5 = m — l — s 

= 2 l ~ 1 ( 2 l - (1 + L)) - (3 l - 2 i-1 (2 + L)) - ( 2 L - (1 + L)) 
= 2 L ~ 1 (2 L - 1) — 3 l + 1 + L, 


using s = dirn(S’) = \Q\ — dim(5 J -) = 2 L — (1 + L). For L = 3,4,5,..., we find 
<5 = 5,44, 259,..., and the deficiency zero and one theorems do not apply. 


More importantly, there exist <5 necessary and sufficient conditions on the 
rate constants for the existence of complex-balancing equilibria, see Horn [9]. 
The conditions involve the Laplacian matrix of the weighted directed graph 
of complexes and reactions, in particular, the quotients of so-called tree 
constants, see 0. For the existence of detailed-balancing equilibria, addi¬ 
tionally the Wegscheider conditions have to be fulfilled, that is, the_products 
of rate constants in a cycle and its reverse must coincide, see ji, 21, 22]. In 
the process of genetic recombination, the rate constants of a reaction and 
its reverse coincide, and all conditions for the existence of complex- and 
detailed-balancing equilibria are fulfilled. 


References 

[1] E. Akin, The Geometry of Population Genetics, vol. 31 of Lect. Notes 
in Biomath., Springer, New York, 1979. 

[2] R. Burger, The Mathematical Theory of Selection, Recombination, 
and Mutation, John Wiley & Sons, 2000. 

[3] A. Dickenstein and M. Perez Millan, How far is complex balanc¬ 
ing from detailed balancing?, Bull. Math. Biol., 73 (2011), pp. 811-828. 






Muller and Hofbauer, 2015 


20 


[4] M. Feinberg, Necessary and sufficient conditions for detailed balanc¬ 
ing in mass action systems of arbitrary complexity, Chemical Engineer¬ 
ing Science, 44 (1989), pp. 1819 - 1827. 

[5] M. Feinberg, The existence and uniqueness of steady states for a class 
of chemical reaction networks, Arch. Rational Mech. Anal., 132 (1995), 
pp. 311-370. 

[6] H. Geiringer, On the probability theory of linkage in Mendelian hered¬ 
ity , Annals Math. Statist., 15 (1944), pp. 25-57. 

[7] A. N. Gorban, General H-theorem and entropies that violate the sec¬ 
ond law, Entropy, (2014). 

[8] J. Higgins, Some remarks on Shear’s Liapunov function for systems of 
chemical reactions, Journal of Theoretical Biology, 21 (1968), pp. 293- 
304. 

[9] F. Horn, Necessary and sufficient conditions for complex balancing in 
chemical kinetics, Arch. Rational Mech. Anal., 49 (1972), pp. 172-186. 

[10] F. Horn and R. Jackson, General mass action kinetics, Arch. Ra¬ 
tion. Mech. Anal., 47 (1972), pp. 81-116. 

[11] L. A. Kun and Y. I. Lyubich, The H-theorem and convergence to 
equilibrium for free multi-locus populations, Kibernetika, (1980), p. 150. 

[12] Y. I. Lyubich, Mathematical Structures in Population Genetics, vol. 22 
of Biomathematics, Springer-Verlag, Berlin, 1992. Translated from the 
1983 Russian original by D. Vulis and A. Karpov. 

[13] L. Markus, Asymptotically autonomous differential systems, in Con¬ 
tributions to the theory of nonlinear oscillations, vol. 3, vol. 36 of Annals 
of Mathematics Studies, Princeton University Press, 1956, pp. 17-29. 

[14] K. Mischaikow, H. Smith, and H. R. Thieme, Asymptotically au¬ 
tonomous semiflows: chain recurrence and Lyapunov functions, Trans. 
Arner. Math. Soc., 347 (1995), pp. 1669-1685. 

[15] S. Muller and G. Regensburger, Generalized mass-action systems 
and positive solutions of polynomial equations with real and symbolic 
exponents (invited talk), in Computer Algebra in Scientific Computing, 
V. P. Gerdt, W. Koepf, W. Seiler, and E. V. Vorozhtsov, eds., vol. 8660 
of Lecture Notes in Computer Science, Springer International Publish¬ 
ing, 2014, pp. 302-323. 

[16] T. Nagylaki, The evolution of multilocus systems under weak selec¬ 
tion, Genetics, 134 (1993), pp. 627-47. 



Muller and Hofbauer, 2015 


21 


[17] T. Nagylaki, J. Hofbauer, and P. Brunovsky, Convergence of 
multilocus systems under weak epistasis or weak selection , Journal of 
Mathematical Biology, 38 (1999), pp. 103-133. 

[18] S. Shahshahani, A new mathematical framework for the study of link¬ 
age and selection , vol. 211 of Memoirs of the AMS, Arner. Math. Soc., 
1979. 

[19] D. Siegel and D. MacLean, Global stability of complex balanced 
mechanisms , J. Math. Chemistry, 27 (2000), pp. 89-110. 

[20] V. M. Vasil’ev, A. I. Vol’pert, and S. I. Hudjaev, The method of 
quasi-stationary concentrations for the equations of chemical kinetics , 
Comput. Math. Math. Phys., 13 (1974), pp. 187-206. 

[21] A. I. Vol’pert and S. I. Hudjaev, Analysis in classes of discontinu¬ 
ous functions and equations of mathematical physics , vol. 8 of Mechan¬ 
ics: Analysis, Martinus Nijhoff Publishers, Dordrecht, 1985. 

[22] R. Wegscheider, Uber simidtane Gleichgewichte und die Beziehungen 
zwischen Thermodynamik und Reactionskinetik homogener Systeme, 
Monatshefte fur Chemie und verwandte Teile anderer Wissenschaften, 
22 (1901), pp. 849-906. 

[23] Y. B. Zel’dovich, The proof of singularity of the solution of mass law 
equations , Zh. Fiz. Khirn., 11 (1938), pp. 685-687. 



